/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
    \file functions.hpp
    \brief Functors Implementation

    This file contains wrappers to perform operetions on \e fun1d and \e fun2d objects.
    \see src/typedefs.hpp
*/

# ifndef FUNCTIONS_HPP
# define FUNCTIONS_HPP 1

# include "typedefs.hpp"
# include "poly.hpp"
# include "point.hpp"

namespace spectral {

/**
    \class TensorProd
    \brief Tensor product

    This will implement the tensor product of two given 1d polynomials:
    \f[ \big(f\otimes g\big)(x,y) = f(x)g(y) \f]
    This class is not meant to be created by the final user directly, but
    is to be used to collect data from the Integrator class.
*/

/**
    \fn TensorProd::TensorProd(Poly const & f,
                   Poly const & g,
                   int idx,
                   int idy,
                   int odx,
                   int ody,
                   double weightX=1.0,
                   double weightY=1.0);
    \brief The constructor
    \param f The x-axis polynomial
    \param g The y-axis polynomial
    \param idx The index of the first polynomial
    \param idy The index of the second polynomial
    \param odx The order of derivation of the first polynomial
    \param ody The order of derivation of the second polynomial
    \param weightX The weight for the first polynomial
    \param weightY The weight for the second polynomial
    We will store a copy of \e f and \e g, because there is an huge
    performance gain with respect to a version of \e TensorProd which
    contains only a reference to \e f and \e g.
*/

/**
    \fn TensorProd::TensorProd(TensorProd const & fg)
    \brief The copy constructor

    This is used to make the TensorProd class compatible with the Product
    class
*/

/**
    \fn TensorProd::TensorProd()
    \brief The default constructor

    This will just call the members default constructors
*/

/**
    \fn TensorProd & TensorProd::operator=(TensorProd const & fg)
    \brief Overloaded assignment operator

    It does what it is supposed to do... make a copy of a TensorProd object
*/

/**
    \fn inline double TensorProd::operator()(double const & x, double const & y) const
    \brief Overloading of the "()" operator

    This is made inline for speed reasons
*/

/**
    \fn Poly TensorProd::operator[](int n) const
    \brief Returns the \e n-sim component

    This will return the first Poly if n=0 else this will return
    the second Poly
*/

/**
    \fn TensorProd _Dx(TensorProd const & fb, double weight)
    \brief Computes the partial derivatives with respect to x

    The \e weight coefficient is used as a moltiplicative coefficient
    for the derivatives. This is used to keep in account the affine
    map between the reference element and the fisical domain and to
    make the derivatives on the fisical domain. This function is
    not intended for the use of the final user, use \e Mesh.Dx instead.
    \see Mesh::Dx
*/

/**
    \fn TensorProd _Dy(TensorProd const & fb, double weight)
    \brief Computes the partial derivatives with respect to y

    The \e weight coefficient is used as a moltiplicative coefficient
    for the derivatives. This is used to keep in account the affine
    map between the reference element and the fisical domain and to
    make the derivatives on the fisical domain. This function is
    not intended for the use of the final user, use \e Mesh.Dy instead.
    \see Mesh::Dy
*/


/**
    \class Functor
    \brief Composition of two functions

    While working with SEM it's always useful to identify a reference
    element with an element inside the computational domain. This is usually
    done defining an affine map: \f$ \phi: \hat{\Omega} \rightarrow \Omega \f$
    and then a given function \f$ f : \Omega \rightarrow R^n \f$ can be composed
    with the map \f$ \phi \f$ so that it can be ``evaluated'' over \f$ \hat{\Omega} \f$.

    This class provides methods and data structures useful to evaluate the composition
    of a given function \f$ f \f$ with a map \f$ \phi \f$. \f$ \phi \f$ must return a
    Point<double> or a double, with an overloaded operator(). If \f$ \phi \f$ returns
    a double then \f$ f \f$ must be a one variable function, otherwise \f$ f \f$ must
    be a two variabile function which returns a double.

    Please notice that there is no easy way to check that the given functors
    have all the required methods implemented, so pay attenction while using this.
*/

/**
    \fn Functor::Functor(T const & f, MAP const & phi)
    \brief The constructor

    T and MAP must be two class with the overloaded operator(). The return
    type of f must be of class Q, while MAP has to return a Point<double>.
*/

/**
    \fn Functor::Functor(Functor<T,MAP> const & func)
    \brief The copy constructor
*/

/**
   \fn inline double Functor::operator()(double const & x) const
   \brief The overloaded operator()

   This will compute:
   \f[ \big( f \circ \phi \big) (x) \f]
*/

/**
    \fn inline double Functor::operator()(double const & x, double const & y) const
    \brief The overloaded operator()

    This will compute:
    \f[ f \big( \phi_1(x,y) , \phi_2(x,y) \big) \f]
*/

/**
    \fn Functor<T,MAP> & Functor<T,MAP>::operator=(
                                Functor<T,MAP> const & func)
    \brief The overloaded operator=

*/


/**
    \class Product
    \brief Product of two function

    This is a wrapper for the product of two functions. The two function
    must have an overloaded operator() with signature:
    \code double operator()(double,double)\endcode
    or:
    \code double operator()(double)\endcode

*/

/**
    \fn Product::Product(T const & t, G const & g)
    \brief The constructor

*/

/**
    \fn Product::Product(Product<T,G> const & fg)
    \brief The copy-constructor

    This is usefull if we want to create an object of the type:
    \code Product<Product<T,G> , G> \endcode to perform an iterated product.
*/

/**
    \fn inline double Product::operator()(double const & x) const
    \brief The overloaded "()" operator

    This is inlined as it's only a wrapper to call the inner functions.
    Please call this only if the functors inside the Product class
    are one variable functions.
*/

/**
    \fn inline double Product::operator()(double const & x, double const & y) const
    \brief The overloaded "()" operator

    This is inlined as it's only a wrapper to call the inner functions.
    Please call this only if the functors inside the Product class
    are two variable functions.
*/

class TensorProd {
    public:
        TensorProd(Poly const & f,
                   Poly const & g,
                   int idx,
                   int idy,
                   int odx,
                   int ody,
                   double weightX=1.0,
                   double weightY=1.0);

        TensorProd();
        TensorProd(TensorProd const & fg);
        TensorProd & operator=(TensorProd const & tp);
        Poly operator[](int n) const;
        inline double operator()(double const & x,double const & y) const;
		inline double operator()(double const & x) const;
        friend class Integrator;
    private:
        Poly _f;             /**< \brief first polynomial */
        Poly _g;             /**< \brief second polynomial */
        int _idx;            /**< \brief index of the first polynomial */
        int _idy;            /**< \brief index of the second polynomial */
        int _odx;            /**< \brief order of the first polynomial */
        int _ody;            /**< \brief order of the second polynomial */
        double _weightX;     /**< \brief weight on the x-axis */
        double _weightY;     /**< \brief weight on the y-axis */
};

template<class T, class MAP>
class Functor {
    public:
        Functor();
        Functor(T const & f, MAP const & phi);
        Functor(Functor<T,MAP> const & func);
        inline double operator()(double const & x) const;
        inline double operator()(double const & x,double const & y) const;
        Functor<T,MAP> & operator=(Functor<T,MAP> const & func);
    private:
        T _f;              /**< \brief Function */
        MAP _phi;          /**< \brief Map */
};

template<class T, class G>
class Product {
    public:
        Product(T const & t, G const & g);
        Product(Product<T,G> const & fg);
        inline double operator()(double const & x) const;
        inline double operator()(double const & x, double const & y) const;
    private:
        T const & _t;              /**< \brief first function */
        G const & _g;              /**< \brief second function */
};

// ====================================================
// End of proper header
// ====================================================

// TensorProd

inline double TensorProd::operator()(double const & x, double const & y) const {
    return _f(x)*_g(y);
}

// Functor

template<class T, class MAP>
Functor<T,MAP>::Functor() {}

template<class T, class MAP>
Functor<T,MAP>::Functor(T const & f,
                        MAP const & phi): _f(f),
                                          _phi(phi) {}

template<class T, class MAP>
Functor<T,MAP>::Functor(Functor<T,MAP> const & func): _f(func._f),
                                                      _phi(func._phi) {}

template<class T, class MAP>
inline double Functor<T,MAP>::operator() (double const & x) const {
    return _f(_phi(x));
}

template<class T, class MAP>
inline double Functor<T,MAP>::operator() (double const & x, double const & y) const {
    Point<double> P(_phi(x,y));
    return _f(P[0],P[1]);
}

template<class T, class MAP>
Functor<T,MAP> & Functor<T,MAP>::operator=(Functor<T,MAP> const & func) {
    _f=func._f;
    _phi=func._phi;
    return *this;
}

// Product

template<class T, class G>
Product<T,G>::Product(T const & t, G const & g): _t(t),
                                                 _g(g) {}

template<class T, class G>
Product<T,G>::Product(Product<T,G> const & fg) : _t(fg._t),
                                                 _g(fg._g) {}

template<class T, class G>
inline double Product<T,G>::operator()(double const & x) const {
    return _t(x)*_g(x);
}

template<class T, class G>
inline double Product<T,G>::operator()(double const & x, double const & y) const {
    return _t(x,y)*_g(x,y);
}

} // namespace

# endif
